import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import os
plt.rcParams['font.family'] = ['SimHei']  # 中文显示
plt.rcParams['axes.unicode_minus'] = False


# 图像读取（彩色和灰度）
def load_images(image_path):
    if not os.path.exists(image_path):
        raise FileNotFoundError(f"图像路径 {image_path} 不存在。")
    color_image = Image.open(image_path).convert('RGB')
    color_array = np.array(color_image)

    gray_image = color_image.convert('L')  # 转为灰度图
    gray_array = np.array(gray_image)

    return color_array, gray_array


# 手动实现二维卷积
def convolution2d_manual(image, kernel):
    image_height, image_width = image.shape
    kernel_height, kernel_width = kernel.shape
    pad_h = kernel_height // 2
    pad_w = kernel_width // 2

    # 对图像进行边界填充，使用镜像填充
    padded_image = np.pad(image, ((pad_h, pad_h), (pad_w, pad_w)), mode='symmetric')

    # 初始化输出数组
    output = np.zeros_like(image, dtype=np.float64)

    # 遍历每一个像素
    for i in range(image_height):
        for j in range(image_width):
            # 提取当前窗口
            window = padded_image[i:i + kernel_height, j:j + kernel_width]
            # 计算卷积
            output[i, j] = np.sum(window * kernel)

    return output


# Sobel算子边缘检测，灰度图
def sobel_edge_detection(gray_image):
    sobel_x = np.array([[-1, 0, 1],
                        [-2, 0, 2],
                        [-1, 0, 1]])
    sobel_y = np.array([[-1, -2, -1],
                        [0, 0, 0],
                        [1, 2, 1]])

    grad_x = convolution2d_manual(gray_image, sobel_x)
    grad_y = convolution2d_manual(gray_image, sobel_y)

    # 计算梯度幅值
    grad_mag = np.sqrt(grad_x ** 2 + grad_y ** 2)

    # 归一化到0-255
    grad_mag = (grad_mag / grad_mag.max()) * 255
    grad_mag = grad_mag.astype(np.uint8)

    return grad_mag


# 给定卷积核的滤波，灰度图像
def given_kernel_filter_gray(gray_image):
    kernel = np.array([[1, 0, -1],
                       [2, 0, -2],
                       [1, 0, -1]])
    convolved = convolution2d_manual(gray_image, kernel)
    # 归一化到0-255
    max_val = np.max(np.abs(convolved))
    if max_val == 0:
        max_val = 1  # 防止除以0
    convolved_norm = (convolved / max_val) * 255
    convolved_norm = np.clip(convolved_norm, 0, 255).astype(np.uint8)
    return convolved_norm


# 颜色直方图
def compute_histogram_color(color_image):
    histograms = []
    for c in range(3):  # 对RGB三个通道
        histogram = np.zeros(256, dtype=int)
        for pixel in color_image[:, :, c].flatten():
            histogram[pixel] += 1
        histograms.append(histogram)
    return histograms


# 纹理特征提取（灰度共生矩阵特征，针对每个彩色通道）
def texture_features_color(color_image):
    gray_level = 256
    texture = []
    for c in range(3):  # 针对RGB三个通道
        glcm = np.zeros((gray_level, gray_level), dtype=int)
        channel = color_image[:, :, c]
        height, width = channel.shape
        for i in range(height):
            for j in range(width - 1):  # 仅考虑水平方向
                pixel_current = channel[i, j]
                pixel_next = channel[i, j + 1]
                glcm[pixel_current, pixel_next] += 1

        # 纹理特征计算
        contrast = 0
        energy = 0
        entropy = 0
        total = np.sum(glcm)
        if total == 0:
            total = 1  # 防止除以0

        for i_val in range(gray_level):
            for j_val in range(gray_level):
                contrast += ((i_val - j_val) ** 2) * glcm[i_val, j_val]
                energy += glcm[i_val, j_val] ** 2
                if glcm[i_val, j_val] > 0:
                    entropy -= glcm[i_val, j_val] * np.log2(glcm[i_val, j_val])

        # 归一化特征
        contrast /= total
        energy /= (total ** 2)
        entropy /= total

        texture.extend([contrast, energy, entropy])
    return np.array(texture)


# 可视化和保存结果
def visualize_and_save(edge_image, filtered_gray_image, histograms, texture, output_path):
    # 显示边缘检测后的灰度图
    plt.figure(figsize=(6, 6))
    plt.imshow(edge_image, cmap='gray')
    plt.title("Sobel算子滤波")
    plt.axis('off')
    plt.show()

    # 保存Sobel算子滤波结果
    edge_image_pil = Image.fromarray(edge_image)
    edge_image_pil.save("sobel_edge.png")
    print("边缘检测结果已保存为 sobel_edge.png")

    # 显示给定卷积核滤波后的灰度图像
    plt.figure(figsize=(6, 6))
    plt.imshow(filtered_gray_image, cmap='gray')
    plt.title("给定卷积核滤波")
    plt.axis('off')
    plt.show()

    # 保存滤波后的灰度图像
    filtered_gray_pil = Image.fromarray(filtered_gray_image)
    filtered_gray_pil.save("filtered_gray_image.png")
    print("滤波后的灰度图像已保存为 filtered_gray_image.png")

    # 颜色直方图
    colors = ['red', 'green', 'blue']
    plt.figure(figsize=(15, 5))
    for i, hist in enumerate(histograms):
        plt.subplot(1, 3, i + 1)
        plt.bar(range(256), hist, color=colors[i].lower())
        plt.title(f"{colors[i]} 直方图")
        plt.xlabel("像素值")
        plt.ylabel("频率")
    plt.tight_layout()
    plt.show()

    # 保存纹理特征
    np.save(output_path, texture)
    print(f"纹理特征已保存至 {output_path}")


# 主函数
def main(image_path, output_path):
    # 加载彩色图像和灰度图像
    color_image, gray_image = load_images(image_path)

    # 1. Sobel算子在灰度图上的边缘检测
    sobel_image = sobel_edge_detection(gray_image)

    # 3. 给定卷积核滤波在灰度图上
    filtered_gray_image = given_kernel_filter_gray(gray_image)

    # 4. 计算颜色直方图
    histograms = compute_histogram_color(color_image)

    # 5. 纹理特征提取
    texture = texture_features_color(color_image)

    # 可视化与保存
    visualize_and_save(sobel_image, filtered_gray_image, histograms, texture, output_path)


if __name__ == "__main__":
    image_path = "图片1.png"
    output_path = "texture_features.npy"
    main(image_path, output_path)
